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The process of stochastic Turing instability on a network is discussed for a specific case study, the 
stochastic Brusselator model. The system is shown to spontaneously differentiate into activator-rich 
and activator-poor nodes, outside the region of parameters classically deputed to the deterministic 
Turing instability. This phenomenon, as revealed by direct stochastic simulations, is explained 
analytically, and eventually traced back to the finite size corrections stemming from the inherent 
graininess of the scrutinized medium. 
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Pattern formation is a rich and fascinating field of in- 
vestigation which extends over distinct realms of applica- 
tions, ideally embracing chemistry, biology, and physics. 
Complex and extremely beautiful patterns can in fact 
spontaneously emerge in spatially extended reaction- 
diffusion systems, as follows a linear instability mecha- 
nism, first described by Alan Turing in a seminal pa- 
per [1| published in 1952. Turing patterns are indeed 
widespread in nature: examples include schemes of auto- 
catalytic reactions with inhibition , the process of bi- 
ological morphogenesis 5 9] and the dynamics of ecosys- 
tems [nUHl • The Turing instability paradigm classically 
relies on mean field, deterministic scenarios. As opposed 
to the usual continuum picture, the intimate discreteness 
of any individual based stochastic models results in finite 
size corrections to the approximated mean-field dynam- 
ics. Under specific conditions, microscopic disturbances 
are enhanced as follows a resonance mechanism and yield 
organized spatio-temporal patterns jl4il71 | . More specif- 
ically, the measured concentration which reflects the dis- 
tribution of the interacting entities (e.g. chemical species, 
biomolecules) can display spatially patched profile, col- 
lective phenomena which testify on a surprising degree 
of macroscopic order, as mediated by the stochastic com- 
ponent of the dynamics. Stochastic Turing patterns [lj|, 
or quasi Turing patterns (l5j |. are found to occur in in- 
dividual based systems, that cannot undergo Turing in- 
stability according to the deterministic reaction-diffusion 
picture. Interestingly, the region of parameters for which 
stochastic patterns are developed is usually larger than 
for conventional Turing patterns, a general observation 
that has been made quantitative for a selection of proto- 
typical case studies. 

Recently, Nakao and Mikhailov [l8| studied the Tur- 
ing patterns formation on large random networks, an im- 
portant direction of investigation presumably relevant in 
e.g. the early stage of the morphogenesis process, since 
morphogens are known to diffuse on the network struc- 
ture of inter-cellular connections. Already in 1971 Oth- 



mer and Scriver [19( investigated the Turing instability 
in network-organized system and developed the needed 
mathematical machineries. Their studies were however 
limited to regular lattice or small networks. By extend- 
ing the analysis to complex heterogeneous network Nakao 
and Mikhailov [l8| opened up the perspective for novel 
applications of the Turing idea to the broad field of the- 
oretical biology [2(J. Applications can be also foreseen 
in other disciplines were network science proves crucial. 
Among others, social studies, in which nodes and links 
are respectively associated to humans and their mutual 
interactions, and the analysis of epidemics spreading, 
which reflects the topological structure of the underly- 
ing mobility networks. 

Starting from this setting, we here propose a gener- 
alization of the work [l^], beyond the deterministic sce- 
nario, by explicitly including the role of demographic, 
finite size fluctuations. In doing so, we will demonstrate 
in this Letter that Stochastic Turing patterns set in, out- 
side the region of parameters deputed to spatial order, as 
predicted within the classical theory based on determin- 
istic reaction-diffusion schemes. 

To this end, we will consider a stochastic version of the 
celebrated Brusselator model which will be placed on 
top of an heterogeneous network of ft nodes. More con- 
cretely, two species, respectively X t and Yi are assigned 
to the generic node i, and therein react according to the 
following chemical reactions (lij : 



A + Ei 
X t + B 
2X t + Yi 



A + X, 
Y t + B 
3X, 

Ei . 



(1) 



The symbol Ei stands for an empty case and formally 
amounts to imposing a finite carrying capacity in each 
node of the network. In other words, we assume that each 
node can host a maximum number N of molecules (or 



2 



agents), including the vacancies [28[. Let us denote by rii 
and nii the total number of elements belonging to species 
X and Y in the node i. The corresponding number of 
empties total hence in N —ni — rrii. The parameters a, b, c 
and d in Eqs. (P) are the reaction rates, while the species 
A and B are enzymatic activators whose concentrations 
are supposed to remain constant during the dynamics. In 
addition to the above activator-inhibitor rules, we assume 
that the molecules can migrate between neighbour nodes 
as dictated by the following reactions: 



with the former, and are given by 



Xi + E 



Ei + X, 



3 7 ^1 ' 

S 



Ei + Y 



(2) 



where /i and S are the diffusion coefficients characteris- 
tic of the two species, and the subscript j denotes the 
generic node connected to i, via the network structure. 
Similar equations governs the diffusion from node j to- 
wards node i. To complete the notation we introduce 
the Sl-dimcnsional vectors n = (m, rii, ...,na) and 
m = (mi, ...,nii, ...,mo) that unequivocally identify the 
state of the system. The process here imagined is in- 
trinsically stochastic. Under the Markov hypothesis, the 
probability P(n, m, t) of seeing the system at time t in 
state (n,m) obeys to a master equation that can be cast 
in the compact form: 



— P(n,m,t) = 53|(«n 4 -l)T(ni + l,mi\m,mi) 

i=l 

+ (e+ - l)T(m - l,mi\m,mi) 

+ {£n, £m, - l)T(m + 1, n%i - l\m,rrn) 

+ (. € ti € mi - l)T(rij - l,rm + \\m,rm) 
n 

+ ^2 Wij [(e^e^ - l)T(m - 1, nj + l\m, rij) 

2=1 

+ (e^.e^i - i)T(nj - l,m + l\ni,nj) 
+ (£,!.,£»,, ~ l)T(mi - 1, rrij + l\rm, rrij) 
+ (d/^ - l)r(mj - 1, rrn + l\rm, rrij)] } 
x P(n,m,t) (3) 



where use has been made of the step opera- 



tors e±/(. 



..,m) 



/(..., Hi ± l,...,m) and 



T(m + l,mi\rii, rrii) 
T(m — 1, mi\rii,mi) 
T(m + 1,171, - l\m, rrn) 
T(m — l,m,+ l\m, rrn) 
T{ni — 1, rij ■ + l\m,rij) 

T(rm — 1, rrij + l\rrn, rrij) 



a N — rii — mi 

n n 
d m 
n n 



n iv 3 
n n 

(i rii N — iij — rtij I 1 



V k + k 3 



S rrii N — rii — ra, / 1 



n n 



N 



where fcj = J2j=i Wij is the degree of the i— th node. 
The factor I /hi + 1/kj takes into account the order of 
selection of species in chemical reactions @. 

The master equation is exact although difficult to 
handle. To progress in the analysis it is customary to re- 
sort to approximated perturbation methods. In the weak 
noise approximation, one can put forward the van Kam- 
pen ansatz [2ll . |22| which, in this context, amounts to 
imposing m/N = (j>i + iu/VN and m, = tpi + &i/VN. 
(j>i and ipi are the mean field concentrations respectively 
associated to the interacting species X and Y. £i, and 
are stochastic fluctuations that originate from finite 
size corrections, normalized by the scalin g fa ctor l/\/~N, 
as dictated by the central limit theorem 21] . For mod- 
erately large system sizes N, the 1/yfN factor is small 
and paves the way to a straightforward perturbative cal- 
culation, generally referred to in the literature as to the 
van Kampen system size expansion. At the leading or- 
der of the perturbative analysis, the mean field equations 
for the deterministic variables are recovered and, for the 
specific problem here investigated, read: 



dr 



n 

3=1 



3 = 1 



dr 



ipi = 9(<l>i,i>i) + 25 



n 

J2 A iJ^3 

3=1 



Q 

^2 A ^3 

3=1 



3=1 



3=1 



(4) 



e„ l ./(n, rrii, ...) = /(n, ...,mj± !,■••), /(•,•) being any where, generalizing the heuristic derivation of [18[, we 



generic function of the state variables. The fl x Q integers 
Wi.j represent the entries of the symmetric adjacency ma- 
trix W, which characterizes the topology of the network. 
Wij is equal to 1 if nodes i and j are connected, and 
otherwise. The transition rates T(n',m'|n, m) link the 
initial state (n,m) to another state (n',m'), compatible 



have introduced the discrete Laplacian A y - = Wij — kiSi. 
with k = Y%=i W l3 and W tJ = (l/h + l/kj)Wij. The 
reaction terms are respectively / = —(& + d)<pi + C(pftpi + 
a(l — 4>i — ipi) and g = btfii — ccpfipi. r is the rescalcd 
time t/(NQ). Cross diffusion terms appear in the ob- 



FIG. 1: The darkened region (yellow on line) in panel (a) 
delineates the Turing instability domain in the (b, c) plane for 
the Brusselator model with o = d = 1, fi = 1 and 8 = 15. 
The (magenta on line) point belongs to the Turing instability 
region and corresponds to b — 76 and c = 950. The (blue 
online) diamond falls outside the region of Turing order and 
is positioned at (76, 1060). In panel (b) the dispersion relation 
© is plotted as a function of both the discrete eigenvalues 
of the network Laplacian (symbols) and their real analogues 
— k 2 (solid line). Circles (magenta online) refer to (b, c) = 
(76,950), while diamonds (bue online) to (6, c) = (76,1060). 
In the analysis we assumed a scale-free network made of 51 = 
200 nodes and mean degree (k) = 20. The network has been 
generated according to the Barabasi- Albert algorithm [25l |. 
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FIG. 2: Simulations of the stochastic chemical model |T}- 
J2} outside the region of Turing order, a = d = 1, b = 76, 
c = 1060, n = 1, 5 = 15. Here N = 1000. The late time 
concentrations per node m/N (resp. rm/N) are plotted in 
the upper panel (resp. lower) panel, as a function of the node 
index i. The (orange online) diamonds are obtained from 
one realization of the stochastic Gillespie algorithm [2(| . The 
network is generated as described in the caption of Fig. [T] 
The stochastic dynamics yields the emergence of two distinct 
activator-rich and activator-poor groups, while the determin- 
istic dynamics is attracted towards the stable (and trivial) ho- 
mogeneous fixed point, dashed (blue online) horizontal line. 



taincd deterministic equations, because of the finite car- 
rying capacity, imposed at the level of the single node 
23[. By relaxing such an assumption 24 1. conventional 
diffusion operators are instead recovered. Similarly, the 
finite carrying capacity assumption reflects in the reac- 
tion contribution a(l — (f>i — tpi) that replaces the usual 
constant term a in the standard Brusselator equations 
[l6|. Although interesting per sc, this modification docs 
not play any substantial role in the forthcoming develop- 
ment: equivalent conclusions can be drawn when working 
in the diluted setting, i.e. away from jamming or crowd- 
ing conditions that inspire the physically sound request 
for a limited capacity to be explicitly accommodated on 
each individual node. 

To look for mean-field Turing instability, one 
needs to introduce a small perturbation to the 
homogeneous equilibrium point (<f>*,i/)*) = ((a + 
^a 2 - 4a6(a + d)/c)/2/(a + d), b/c/<j>*) of the determin- 
istic system Q and carry out a linear stability analysis. 
In formulae, (^j,Vi) = (</>* + Stpi). Following 

fl8j , and to exploit the linearity of the resulting equations 
for the perturbation amounts [2j|, we find it convenient 



to expand Stfii and 6(j>i as: 

Q 

,AaT„,(") 



(5) 



Q = l 



where v^ Q ' = (v[ a \ . . . , v^) stand for the eigenvectors 
of the Laplacian operator corresponding to the eigenvalue 
A« 0. 

By inserting Eqs. ([5]) into the linearized differential 
equation for the perturbations 5<pi and Sipi, one obtains 
the usual characteristic equation for A Q , which can be 
here cast in the form: 



dct 



+ <5i/j*A q 



g^ + S(l — 4>*)^a - A Q 







(G) 



where f q = df /dq and g q = dg/dq for q = <j>, rp. 

The Turing instability occurs, and the perturbation 
gets thus amplified, if A a (A a ) is positive for some value 
of A Q . In this respect, and as already remarked in 
A a plays the role of —k 2 for continuous media, where 
k stands for the wavenumber of the plane wave mode. 
In Fig. [ljb), the dispersion relation is plotted for two 
distinct choices of the parameters (see legend) . Symbols 



4 



refer to the discrete linear growth rates X a , as function of 
the corresponding Laplacian eigenvalues A Q . The solid 
line represents the homologous dispersion relations, as 
obtained working within the continuous representation 
(A Q — > — k 2 ). The upper curve (panel (b) of Fig. [TJ 
circles) signals the presence of an instability. A signifi- 
cant fraction of the discrete rates Aq, is in fact positive. 
Conversely, the other profile (diamonds) is obtained for 
a choice of the chemical parameters that yields linear 
stability. By tuning the parameters, and evaluating the 
corresponding dispersion relation, one can eventually sin- 
gle out in a reference parameter space the region deputed 
to the instability. This is done in Fig. \T[&) working in 
the plan (6, c): the region of Turing instability, as pre- 
dicted by the deterministic analysis, is filled with a uni- 
form colour (yellow online). The (blue online) diamond 
falls outside the region of Turing order and points to the 
parameters employed in depicting the stable dispersion 
curve in panel (b). Similarly, the circle (magenta on- 
line) refers to the unstable profile. In this latter case, 
performing a direct integration of the mean-field equa- 
tions ((3]) one observes the spontaneous differentiation in 
activator-rich and activator-poor groups, as discussed in 
[l8j]. A stochastic simulation can be also carried out, 
using an ad hoc implementation of the Gillespie Monte 
Carlo scheme [26|. In the movies annexed as supple- 
mentary material, the two dynamics, deterministic vs. 
stochastic, are compared. Finite size fluctuations materi- 
alize in a modest perturbation (oc 1/y/N) of the idealized 
mean-field dynamics. 

Substantially different, is instead the scenario that 
is eventually recovered when comparing the simulations 
outside the region deputed to Turing instability. Setting 
the parameters (b — 76, c = 1076) to the values (b = 76, 
c = 1060) that correspond to the diamond (blue online) 
of Fig. [Ha), the deterministic simulations always con- 
verge to the homogeneous fixed point, the concentrations 
of the species being therefore identical on each node of 
the network. At variance, a fragmentation into distinct 
groups is clearly observed in the stochastic simulations. 
The late time evolution of the stochastic system, as com- 
pared to the corresponding (trivial) deterministic solu- 
tion, is displayed in Fig. [2 The effect of the stochastic 
driven polarization can be further realized when inspect- 
ing the annexed movies, which enables one to appreciate 
the full time evolution of the discrete dynamics. As for 
the case of continuous media, the endogenous stochastic 
noise is amplified and drives the formation of spatially 
extended, self-organized patterns outside the region of 
classical Turing order. Following [la ], we call these self- 
organized, asymptotically stable configurations, Stochas- 
tic Turing patterns on a network. 

To gain analytic insight into the above mechanism, one 
can return to the van Kampen perturbative analysis and 
consider the next to leading approximation. One obtains 
a system of Langevin equations [l4, 16 1 for the fluctua- 




FIG. 3: Power Spectrum of fluctuations for species X as a 
function of A a , lj = (symbols). The solid line is the power 
spectrum calculated for a continuum media, i.e. when A a 
is replaced by — k 2 , where k denotes the wavenumber of the 
plane wave mode. The curve refers to a = d = 1, b = 76, c = 
1060, /it = 1, 5 = 15, a choice of parameters that correspond to 
operate outside the region of Turing instability (diamond in 
[TJa)). The network is constructed as specified in the caption 
of Fig. ED 



tions £ S i, s = 1, 2: 
dr 



where r) S i is a Gaussian noise with zero mean and corre- 
lator given by (r]s,i(T)ri r j(T')) = B sr ^j5 TT i. The explicit 
form of the matrices M. and B will be given in the Ap- 
pendix. Define then following transformation: 



. e -J-T>0 



(8) 



j denoting here the imaginary unit. The above operation 
is inspired to the Fourier transform: instead of expand- 
ing on a basis of plane waves, it is here natural to project 
the fluctuations along the fl independent directions rep- 
resented by the eigenvectors of the discrete network 
Laplacian. One can hence define a power spectrum of 
fluctuations of species s = 1,2, P s (uj,A a ) = (|£ s | 2 ), in 
completely analogy with what it is customarily done in 
conventional Fourier analysis. In practical terms, the 
generalized power spectrum P s (uj, A q ) quantifies the por- 
tion of the signal power that is associated to given time 
(a;) or spatial frequencies (A Q ) range. Some details of the 
calculations are given in the Appendix. 

In Fig. [3] the analytical power spectrum of species X, 
is plotted in the plane u = 0, as a function of A Q , for the 
parameters selection that corresponds to the diamond in 
Fig. ([T]), i.e. outside the region of deterministic Tur- 
ing instability. Symbols are obtained by sampling the 
power spectrum over the discrete Laplacian eigenvalues 
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A a . The solid line stands for the power spectrum cal- 
culated in the continuous limit when the discrete A Q is 
replaced with its real counterpart — > — k 2 (see Appendix 
for a discussion related to this point). A clear peak is dis- 
played [3l| , a finding that explains in turn the outcome of 
the stochastic based simulations reported in Fig. [5J prov- 
ing on a formal ground that stochastic Turing patterns 
do exist on a network topology. 

In conclusion we have considered in this Letter the 
stochastic dynamics of the Brusselator model on a net- 
work. The model, representative of a broad class of sys- 
tems that display Turing order in the mean field limit, 
has been investigated both analytically and numerically. 
In particular, we provide evidences on the intrinsic abil- 
ity of the system to develop spatially heterogeneous con- 
figurations, outside the region of parameters classically 
deputed to Turing (deterministic) order. These self- 
organized patterns, reminiscent of the Turing instability, 
result from the spontaneous amplification of the demo- 
graphic noise, stemming from the intimate discreteness of 
the scrutinized medium. Our analysis extends therefore 
the concept of Stochastic Turing order to the vast domain 
of network science, a discipline of paramount importance 
and cross-disciplinary interests. Further investigations 
can be planned working along these lines, and so explore 
the surprising degree of macroscopic order that can even- 
tually originate from the noisy microscopic dynamics, for 
stochastic based systems defined on a network topology. 
As an example, stochastic travelling waves can be imag- 
ined to occur as a natural extension of [24[ . Incidentally, 
we also emphasize that the finite carrying capacity mech- 
anism here imposed at the node level could turn useful 
to model those phenomena where jamming on a network 
topology are expected occur. 



ing entries for matrix M. ^ NS ^ : 

M { ^ S) = -a-b-d + 2c<j>*^*, 
M 12 = -a + c<p , 
M ( 21 S) = b-2c<j>*ip*, 

. A (NS) ,* 2 

M\ 2 = -af> . 
The elements of matrix A4^ SP > read instead: 
M^ P) = 2 M (1-V*), 

m[ s 2 p) = 2n4r, 

M { i P) = 26p, 

M ( i p) = 26(1 -<n 

As concerns the matrix B. one finds: 



B\\,a 


= D 1 + 


Bl2,ii 


= £>21,m 




= D 2 + 


Bii,« 


= ~ {h 




= Bai,« 




= Vfc 



Hfl2 



where: 



Dt = a(l-<f>* -V>*) + <£* {b + ccj>*4>* +d) 

Fi = 4/z<£*(l - </>* -V*) 

C = -<t>* (6 + cW) 

D 2 = 4>* (b + c^r) 

H 2 = 45ip*(l -(/)*- ip*) 



At the next-to-leading order approximation in the van 
Kampcn expansion, one obtains a Fokkcr Planck equa- 
tion for the probability distribution of fluctuations, which 
is equivalent to the Langevin equation (|7|). For each 
fixed pair of nodes, i, j, the 2x2 matrix A^ sr ,ij may 
be decomposed as the sum of two contributions, one 
relative to the activator-inhibitor reactions (non spa- 
tial components (NS)), and the other associated to the 
diffusion process (spatial components (SP)): Ai sr .ij = 
M sT ( N S) _|_ _A4 sr ( SP ) The above matrices are eval- 
uated at the mean-field fixed points </>*, tp* ■ 

After a lengthy calculation [l6| one obtains the follow- 



Matrix B can be hence cast in the equivalent form: 
B SS)ij = (D s + kiH,)Sij - Q- + WijH s 

Brs,ij — B sr ^[j — CSij (9) 

for r, s = 1, 2. 

Performing the transformation (|5J| on both sides of the 
Langevin equation (|7|), we get 

2 
r=l 

where the term fjf denotes the transform of the noise. 
Introducing the matrix 

<& sr = jL05 sr - (Mg 5) + X^ P) Aa) (11) 
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we get = X)r=i ®srV? > an( i thus the power spectrum 
of the s-th species is given by 

2 

P s (u,A a ) = (|C| 2 ) = J2 *^<w)(*L) _1 ( 12 ) 



\fc=l 



It can be shown that (77™ t)£) = ^ . B r fc,ijuj a ' 1^ . An 
explicit form for the 2x2 matrix (fj"fj^) can be derived 
by making use of Eqs Let us focus on the non trivial 
contribution r = k: 

/ , °rr,ij <Aj — 



hJ 



E 



Hr A h 



(a) (a) 



(a) (a) 

V; 



,(«)„,(«) 



J u i 



where in the last step we made use of J^. v\ v\ = 1 and 
y~]j AijVj = A a vf. The component r 7^ k is trivially 
equal to C. 

The power spectrum is fully specified as function of A a 
and ui. Figure [3] is obtained by setting uj = in the above 
formulae. We emphasize that the same result can be 
recovered from the continuum medium power spectrum 
[l6l | provided — k 2 is replaced by the discrete eigenvalue 
A Q . 

We end this Appendix by providing a list of explana- 
tory captions to the movies annexed as supplementary 
material. 

1. mf_inside_X.mov and mf-inside-Y.mov show re- 
spectively the time evolution of the mean field con- 
centrations 4>i and ipi- The data are obtained by 
integrating the governing mean field equations (@| . 
Parameters are chosen so to yield a Turing insta- 
bility (magenta circle in Fig. [1]). At time t = 
a small perturbation is applied to perturb the ho- 
mogeneous fixed point. Then, the system evolves 
toward a stable, non-homogeneous stationary con- 
figuration. 

2. st-inside-X.mov and st-inside-Y.mov show the re- 
sult of the stochastic simulations (fluctuating blue 
circles), for respectively n.i/N and rrii/N. The red 
symbols refer to a late time snapshot of the deter- 
ministic dynamics. Parameters are set as specified 
above: the system is hence inside the region of Tur- 
ing instability. Notice that the noise that takes the 
system away from the homogeneous fixed point is 
now endogenous to the system and not externally 
imposed. 

3. st-outside-X.mov and st-outside-Y.mov report on 
the stochastic simulations (blue symbols), for re- 
spectively rii/N and rrii/N. The parameters are 



now assigned so to fall outside the region of Turing 
order (blue diamond in Fig. [T]). The mean field 
solutions are not destabilized and converge to the 
stable fixed point. At variance, the stochastic sim- 
ulations evolve toward a non homogeneous state. 
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the code to generate the network. 
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